library(MASS)
library(gplots)
library(ggplot2)
library(Hmisc)
library(sciplot)
library(RItools)


rm(list=ls())

qtdata <- read.csv(file="Rrep_2013_07_29.csv",head=TRUE,sep=",")
qtdata2 <- qtdata[ which(qtdata$treatment2!="NA" & qtdata$approve_dum !="NA"), ]
attach(qtdata2)

#meanbyti<- as.matrix(aggregate( approve_dum ~ treatment2, data=qtdata2, FUN=mean, na.rm=TRUE))
nbyti<- as.matrix(aggregate(approve_dum ~ treatment2, data=qtdata2, FUN=function(x) sum(!is.na(x)) ))

n0<-nbyti[1,2]
n1<-nbyti[2,2]
n2<-nbyti[3,2]
n3<-nbyti[4,2]

qtdata2_app <- qtdata2[ which(qtdata2$approve_dum =="1"), ]

abyti<- as.matrix(aggregate(approve_dum ~ treatment2, data=qtdata2_app, FUN=function(x) sum(!is.na(x)) ))


a0<-abyti[1,2]
a1<-abyti[2,2]
a2<-abyti[3,2]
a3<-abyti[4,2]

n<-5000

pi_0 <- as.matrix(rbeta(n, a0+0.5, n0-a0+0.5))
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))

t0 <- as.matrix(rep(0,n, nrow=n, ncol=1))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))

tis <- as.matrix(rbind(t0,t1,t2,t3))
pis <- as.matrix(rbind(pi_0,pi_1,pi_2,pi_3))

betas <- as.data.frame(cbind(tis,pis))
betas$treatment2name = factor(betas$V1, labels = c("Null","Int. Agr.","Placebo","Econ."))


###
# Figure 1
###
#This is the main plot, approval ratings by treatment group

lineplot.CI(x.factor = betas$treatment2name, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "Approval %",
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )

###
# Set Up for Figures 2-3
###
# Treatment Effects, broken down by trade preferences
###

### IL Beta Distributions ###
qtdata2_il <- qtdata2[ which(qtdata2$treatment2=="2"), ]

# Number of respondents for IL treatment, by trade category
nbypi_il<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_il, FUN=function(x) sum(!is.na(x)) ))

np1_il<-nbypi_il[1,2]
np2_il<-nbypi_il[2,2]
np3_il<-nbypi_il[3,2]

# Number of approvals for IL treatment, by trade category
qtdata2_il_app <- qtdata2_il[ which(qtdata2_il$approve_dum=="1"), ]

abypi_il<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_il_app, FUN=function(x) sum(!is.na(x)) ))

ap1_il<-abypi_il[1,2]
ap2_il<-abypi_il[2,2]
ap3_il<-abypi_il[3,2]

# Beta draws with data generated parameters
pi_1_il <- as.matrix(rbeta(n, ap1_il+0.5, np1_il-ap1_il+0.5))
pi_2_il <- as.matrix(rbeta(n, ap2_il+0.5, np2_il-ap2_il+0.5))
pi_3_il <- as.matrix(rbeta(n, ap3_il+0.5, np3_il-ap3_il+0.5))

# Indicators for trade category
p1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
p2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
p3 <- as.matrix(rep(3,n, nrow=n, ncol=1))

# Indicators for IL treatment
t2b <- as.matrix(rep(2,n*3, nrow=n, ncol=1))

ps <- as.matrix(rbind(p1,p2,p3))
pis_il <- as.matrix(rbind(pi_1_il, pi_2_il, pi_3_il))

# Matrix of betas for IL treatment
betas_il <- as.data.frame(cbind(t2b,ps, pis_il))
betas_il$tradecategname = factor(betas_il$V2, labels = c("Pro Free Trade","No Opinion","Protectionist"))



### NL Beta Distributions ###
qtdata2_nl <- qtdata2[ which(qtdata2$treatment2=="1"), ]

# Number of respondents for NL treatment, by trade category
nbypi_nl<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_nl, FUN=function(x) sum(!is.na(x)) ))

np1_nl<-nbypi_nl[1,2]
np2_nl<-nbypi_nl[2,2]
np3_nl<-nbypi_nl[3,2]

# Number of approvals for NL treatment, by trade category
qtdata2_nl_app <- qtdata2_nl[ which(qtdata2_nl$approve_dum=="1"), ]

abypi_nl<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_nl_app, FUN=function(x) sum(!is.na(x)) ))

ap1_nl<-abypi_nl[1,2]
ap2_nl<-abypi_nl[2,2]
ap3_nl<-abypi_nl[3,2]

# Beta draws with data generated parameters
pi_1_nl <- as.matrix(rbeta(n, ap1_nl+0.5, np1_nl-ap1_nl+0.5))
pi_2_nl <- as.matrix(rbeta(n, ap2_nl+0.5, np2_nl-ap2_nl+0.5))
pi_3_nl <- as.matrix(rbeta(n, ap3_nl+0.5, np3_nl-ap3_nl+0.5))

# Indicators for trade category
# [Don't need to repeat]

# Indicators for NL treatment
t1b <- as.matrix(rep(1,n*3, nrow=n, ncol=1))

# ps <- as.matrix(rbind(p1,p2,p3))  #[don't need to repeat]
pis_nl <- as.matrix(rbind(pi_1_nl, pi_2_nl, pi_3_nl))

# Matrix of betas for NL treatment
betas_nl <- as.data.frame(cbind(t1b,ps, pis_nl))
betas_nl$tradecategname = factor(betas_nl$V2, labels = c("Pro Free Trade","No Opinion","Protectionist"))




### EC Beta Distributions ###
qtdata2_ec <- qtdata2[ which(qtdata2$treatment2=="4"), ]

# Number of respondents for EC treatment, by trade category
nbypi_ec<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_ec, FUN=function(x) sum(!is.na(x)) ))

np1_ec<-nbypi_ec[1,2]
np2_ec<-nbypi_ec[2,2]
np3_ec<-nbypi_ec[3,2]

# Number of approvals for EC treatment, by trade category
qtdata2_ec_app <- qtdata2_ec[ which(qtdata2_ec$approve_dum=="1"), ]

abypi_ec<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_ec_app, FUN=function(x) sum(!is.na(x)) ))

ap1_ec<-abypi_ec[1,2]
ap2_ec<-abypi_ec[2,2]
ap3_ec<-abypi_ec[3,2]

# Beta draws with data generated parameters
pi_1_ec <- as.matrix(rbeta(n, ap1_ec+0.5, np1_ec-ap1_ec+0.5))
pi_2_ec <- as.matrix(rbeta(n, ap2_ec+0.5, np2_ec-ap2_ec+0.5))
pi_3_ec <- as.matrix(rbeta(n, ap3_ec+0.5, np3_ec-ap3_ec+0.5))

# Indicators for trade category
# [Don't need to repeat]

# Indicators for EC treatment
t4b <- as.matrix(rep(4,n*3, nrow=n, ncol=1))

# ps <- as.matrix(rbind(p1,p2,p3))  #[don't need to repeat]
pis_ec <- as.matrix(rbind(pi_1_ec, pi_2_ec, pi_3_ec))

# Matrix of betas for EC treatment
betas_ec <- as.data.frame(cbind(t4b,ps, pis_ec))
betas_ec$tradecategname = factor(betas_ec$V2, labels = c("Pro Free Trade","No Opinion","Protectionist"))


### PL Beta Distributions ###
qtdata2_pl <- qtdata2[ which(qtdata2$treatment2=="3"), ]

# Number of respondents for PL treatment, by trade category
nbypi_pl<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_pl, FUN=function(x) sum(!is.na(x)) ))

np1_pl<-nbypi_pl[1,2]
np2_pl<-nbypi_pl[2,2]
np3_pl<-nbypi_pl[3,2]

# Number of approvals for PL treatment, by trade category
qtdata2_pl_app <- qtdata2_pl[ which(qtdata2_pl$approve_dum=="1"), ]

abypi_pl<- as.matrix(aggregate(approve_dum ~ tradecateg, data=qtdata2_pl_app, FUN=function(x) sum(!is.na(x)) ))

ap1_pl<-abypi_pl[1,2]
ap2_pl<-abypi_pl[2,2]
ap3_pl<-abypi_pl[3,2]

# Beta draws with data generated parameters
pi_1_pl <- as.matrix(rbeta(n, ap1_pl+0.5, np1_pl-ap1_pl+0.5))
pi_2_pl <- as.matrix(rbeta(n, ap2_pl+0.5, np2_pl-ap2_pl+0.5))
pi_3_pl <- as.matrix(rbeta(n, ap3_pl+0.5, np3_pl-ap3_pl+0.5))

# Indicators for trade category
# [Don't need to repeat]

# Indicators for PL treatment
t3b <- as.matrix(rep(3,n*3, nrow=n, ncol=1))

# ps <- as.matrix(rbind(p1,p2,p3))  #[don't need to repeat]
pis_pl <- as.matrix(rbind(pi_1_pl, pi_2_pl, pi_3_pl))

# Matrix of betas for PL treatment
betas_pl <- as.data.frame(cbind(t3b,ps, pis_pl))
betas_pl$tradecategname = factor(betas_pl$V2, labels = c("Pro Free Trade","No Opinion","Protectionist"))


###
# IL vs NL, by trade preferences
###

# Making large betas matrix
betas_ilnl <- as.data.frame(rbind(betas_il, betas_nl))

# Approval by tradecateg, grouped by treatment (IL NL)
betas_ilnl$tmnames = factor(betas_ilnl$V1, labels = c("Null","Int. Agr."))

###
# PL vs NL, by trade preferences
###

# Making large betas matrix
betas_plnl <- as.data.frame(rbind(betas_pl, betas_nl))

# Approval by tradecateg, grouped by treatment (PL NL)
betas_plnl$tmnames = factor(betas_plnl$V1, labels = c("Null","Placebo"))

# Making large betas matrix
betas_ecnl <- as.data.frame(rbind(betas_ec, betas_nl))

# Approval by tradecateg, grouped by treatment (EC NL)
betas_ecnl$tmnames = factor(betas_ecnl$V1, labels = c("Null","Econ."))


###################
# Combo Figure (3)
###################
par(mfrow=c(2,1),oma=c(0,0,0,0), mar=c(2.9,4.1,0.9,0.9))

lineplot.CI(x.factor = betas_plnl$tradecategname, response = betas_plnl$V3, group=betas_plnl$tmnames, type="p",
	data = betas_plnl, main = "", xlab = "Trade Preferences", ylab = "Approval %",
	col = c("red","blue"),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)),
	x.leg=2.7, y.leg=0.59, ylim=c(0.40,1))

lineplot.CI(x.factor = betas_ecnl$tradecategname, response = betas_ecnl$V3, group=betas_ecnl$tmnames, type="p",
	data = betas_ecnl, main = "", xlab = "Trade Preferences", ylab = "Approval %",
	col = c("red","blue"),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)),
	x.leg=2.7, y.leg=0.59, ylim=c(0.40,1))

dev.print(postscript, "C:/WTO_Survey/Drafts/approve_by_ft_combo.eps", horizontal=FALSE)




##################
# R1: Bowers/Hansen Covariate Balance
##################
xBalance(tmt_ilaw~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum, data=qtdata2,
	report="all")
xBalance(tmt_econ~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum, data=qtdata2,
	report="all")
xBalance(tmt_plac~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum, data=qtdata2,
	report="all")
xBalance(tmt_null~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum, data=qtdata2,
	report="all")


### Full barrage of covariates
xBalance(tmt_ilaw~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum + pk_sum + iso_standard + ethno_standard + working + ab_med_income + repub_dum + conserv_dum + union_dum, data=qtdata2,
	report="all")
xBalance(tmt_econ~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum + pk_sum + iso_standard + ethno_standard + working + ab_med_income + repub_dum + conserv_dum + union_dum, data=qtdata2,
	report="all")
xBalance(tmt_plac~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum + pk_sum + iso_standard + ethno_standard + working + ab_med_income + repub_dum + conserv_dum + union_dum, data=qtdata2,
	report="all")
xBalance(tmt_null~ age + male + race_white + race_black + race_hisp + race_asian + married_dum + college_dum + pk_sum + iso_standard + ethno_standard + working + ab_med_income + repub_dum + conserv_dum + union_dum, data=qtdata2,
	report="all")

